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I. Introduction 


Transport, along with such processes as radiation and wall recycling, 
determine the energy and particle confinement times as well as density 
and temperature profiles for ions and electrons. As a result, trans- 
port and atomic processes determine whether a particular equilibrium 
is accessible and how it will evolve. Transport predicts how energy 
and particle confinement times scale with average density, gross 
energy, radius, and temperature of the plasma. These so-called scaling 
laws are an important guide for the design of new confinement devices. 

Transport codes are one of the essential research tools with which 
experimental data can be interpreted and theoretical ideas can be 
quantitatively tested. They isolate the role of various processes 
involved in transport. As an illustration, on the next page Fig, 1 
shows the relative importance of the anomalous electron heat and plasma 
transport when the neoclassical heat transport is switched off. At 
the present time transport codes are much more useful as research tools 
than as predictive, reliable engineering design tools. This is because 
the transport rates are such as not to make the codes predictive. As 
a result, many of the major codes are not of completely general utility 
but have been developed with special emphasis on particular in-house 
experimental needs. This is no reflection on the effectiveness of the 
codes since simple treatments are quite often adequate. An example of 
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a research code which was designed, among other things, to test the 
range of validity of Pfirsch-Schl uter diffusion is given in Appendix III. 
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FIG.. 1* Density dependence of total energy confinement time for 
deuterium discharges in ISX-A with X « 120 kA and B t « 
12.8 kG. The solid line shows full transport model pre- 
dictions, fit to the experimental point at n *» 1.5 * 10 13 
cm • The- dashed line shows the corresponding results 
with ion neoclassical thermal conductivity neglected in 
the transport model. 


* 

R* E, Walt 2 and G, E, Guest, GA-A15048, to be published. 
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II, Summary 

A workshop on 1-D and 2-D transport in tokamaks was held at Ithaca, 

New York on August 2-4, 1978. The purpose of the Workshop was to 
assess the status of physical models used in transport calculations 
and to evaluate the maturity of 2-D transport codes in predicting 
operating parameters of such confinement devices as Alcator, PLT, 

Doublet III, and TFTR. 

The Workshop consisted of in-depth discussions of the following 
topics: 

Status of 1-D codes, problems where 2-D treatment is necessary or use- 
ful, status of the treatment of fundamental processes, successful 
models, boundary and wall effects, 3-D and velocity space effects, and 
numerical algorithms used in transport codes. 

Below is a brief summary of some of the discussions and conclusions: 

Transport Modules . Since transport of ionized species along field 
lines is much faster than across field lines, transport in two and even 
three dimensions can be reduced to a one-dimensional problem dealing 
with transport across flux surfaces. Both 1-D and 2-D, or so-called 
l^a-D transport modules calculate radial, i.e., one-dimensional transport. 
The difference between the two classes of codes lies in how they treat 
the coupling between flux surface motion and conservation of mass, 
momentum, and flux. 



4 


1- D Codes . The surfaces of constant poloidal flux are assumed to be 
concentric toroids of circular cross section. The fact that experimental 
profiles are roughly symmetric justifies this assumption. It is 
important to emphasize that existing 1-D codes may contain a great deal 
of 2-D and 3-D information. For example, neutral beam deposition is 
calculated on the basis of 2-D particle orbits, fast ions are calculated 
on the basis of three velocity variables, flux surface averages are 
based on 2-D flux surfaces. All of this information is folded as 
coefficients into a differential equation that depends on radius and 
time only. 

2- D Codes . Even though the transport module is 1-D, these codes 
solve fully 2-D problems to determine magnetic flux surfaces. Flux 
surface averages are performed on 2-D flux surfaces whose shape can 
change as a result of transport. 

2-D codes must be used when there are moderate changes in magnetic 
geometry, for the design of magnetic circuitry, in MHD stability 
studies and ion beam deposition studies. It must be emphasized that 
the efficient use of 2-D codes increases CPU time by only 50%. 

Physics Models . It was agreed that the models dealing with atomic 
processes and neutral particles are well in hand. Non-coronal models 
arp ava-iiahio are not yet in general use. 


beam deposition, alpha particle heating, plasma 
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compression, and general gross 2-D motion of magnetic flux surfaces 
have reached advanced stages of development. 

The physical boundary region of the plasma is not understood quanti- 
tatively. This includes production and inward transport of recycled 
material, partition of energy flow between wall and limiter, and non- 
coronal radiation from impurities. 

There is little quantitative understanding of the anomalous electron 
heat and plasma transport. The simplest phenomenological approach 
has been to assign a value to the anomalous heat conductivity inversely 
proportional to the density. This gives the Alcator scaling, but 
the variation of conductivity by a factor of 3 to 5 from machine to 
machine has not been correlated with any systematic dependence. 

More work needs to be done on anisotropic and other types of distortion 
of plasma equilibrium by neutral beam injection. 

Future Directions 

The development of transport codes will be affected by the amount and 
quality of experimental data available. 

In the near future the codes will undergo improvements that will include 
the coupling of large impurity diffusion calculations to transport codes, 
and development of 3-D neutrals models to calculate effects of gas 
puffing. The coupling of MHD stability codes with transport calculations 
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will give information on the accessibility of high-beta states and on 
the effect of ballooning modes in these states. This area of research 
will receive increased emphasis at each major laboratory. 

As the use of more physics and of higher spatial and velocity dimensions 
increases in the future, the codes will become more specialized to 
particular problems or machines. At the same time, development of 
large complex "master" codes will continue. These will be primarily 
used to check the output of the more specialized codes. 


III. Recommendations 


The Workshop participants have identified two areas that need particular 
attention. These are the lack of an experimental tokamak data base 
and the unavailability of neoclassical transport coefficients for 
arbitrary coll isional ity, geometry, and aspect ratio. 

A. It is recommended that a Committee on the Tokamak Data Base 
be established to develop procedures by which properly aged data on 
tokamak experiments be made available to code builders. The Committee 
should consist predominantly of experimentalists, with a participation 
of theoreticians working in the field of transport. 

B. It is recommended that a compendium of transport coefficients 
be made available in the parameter regions where they are understood. 

It is also recommended that encouragement be given to develop formulas 
for coefficients in parameter regions where they are not so well under- 
stood. For example, work should be encouraged on such things as neo- 
classical ion thermal conductivity and trapped particle corrections 


to the Ohm's law resistivity. 
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I , F hy s i cal Mod els in Transport Codes 

hi , rdcr to calculate transport in a fusion device, a code must contain 
dealing with transport of several plasma species, anomalous 
•'Treats, the effect of impurities, radiation losses, neutral injection, 
ftdupnic compression, alpha particle heating, and plasma-wall and plasma- 
t no ter Interaction, The code must be such as to ensure particle and 
n-rgy balance. Below we will discuss briefly each of these effects. 


1 i:i ,,, 3 .Jpecies. In order to know profiles in a plasma, it is necessary 
to Know how each species diffuses across the magnetic field. Consequently, 
transport codes treat wore than one species of particles and more than 
one charge state of each species. 


C^llousjvansport. In addition to classical transport due to Coulomb 
collisions and electric and magnetic field ^homogeneities, a plasma 

" M/ aC, '° SS magnetic F1eld to various types of collective 

instabilities that are largely unrelated to classical transport. Such 

transport is called anomalous and for electrons it is in general higher 

cnan classical transport. A discussion of U i 

scussion of the problems in calculating 

anomalous transport is given in Section VI. 

InountxJifects and RadUti 

- Impurities affect tokamaks by 

eigy, contributing to z gff and thus increasing the ohmic 

^ Snd ^ dl,Ut1n9 ^ P,aSma ' The calculation of impurity 
* 9reat Pr ° 9ress ln the last few years. Atomic rates 
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for computing the impurity radiation power losses have been calculated 

for essentially all the known elements. The most common approximation 

used is to assume "coronal equilibrium" (ionization and recombination 

times short compared to diffusion times). This approximation is 

reasonably good for high Z elements (Z > 70), especially near the center 

of a tokamak. It is convenient since the power radiated can be computed 

from a simple formula P 7 = n n L,(Te) where P_ is the power radiated in 
3 

watts/cm , n g is the electron density, n 2 is the impurity density, and 
L 2 (Te) is the luminosity of the element which is a function of the electron 
temperature. Tables of fits to L z (Te) have been compiled and are used in 
most transport codes. 

Coronal equilibrium is not a good approximation for low Z elements 

in devices with T o > 200 eV. In these machines, the low Z impurities 

recycle at the edge. For example, oxygen ions enter as nearly neutral 

0 , etc., and are ionized rapidly to 0 , and eventually to 0 and 
+8 

0 . The partially ionized ions radiate furiously while they are being 

ionized. A single oxygen ion may radiate as much as 10 keV while being 
ionized (compared to ^ 400 eV for the sum of ionization potentials for 

in if 0.0 

0 to 0 ). Once ionized the 0 migrates into the center and comes out 

again, persisting almost to the limiter or wall, where it recombines and 
recycles again. 


Models for this radiation have been constructed and used in some trans- 
port codes, but are not in common usage today. 
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Neutra l Injection , As the beam particles penetrate the plasma, they 
undergo charge exchange with plasma ions and become ionized. In order 
to calculate the deposition of energy by neutral beam during injection 
beam ion trajectories must be obtained. These are obtained either by 
numerical integration of the standard guiding center equations for a 
torus or by calculating the effect of Coulomb collisions through the 
use of the Fokker-Planck equation. Neutral injection codes must be 
able to calculate deposition not only by full energy but also by 1/2 
and 1/3 energy components of the beam. 


Neutral .Tr ansport . Since particle confinement times in tokamaks are 
much shorter than the average pulse length, particle sources are 
important. The model used in calculations assumes that plasma diffuses 
out to the limiter and is neutralized. It recycles into the discharge 
as neutral gas which is ionized to form density sources. The neutral 
gas atoms can charge exchange with the plasma ions, and thus penetrate 
quite far into the plasma. Other sources of neutral gas beside recycl Ing 
Plasma are gas puffing, recombination, and neutral beam deposition. 


Neutral gas transport is calculated from the linear Boltzmann equation 
Semi -analytic slab models, discrete ordinate techniques, and Monte 
Carlo methods have been employed in the neutral gas problem. 


Mimioc^res^. The mgnetic m m Qf a cha ^ d part . cie 
is invariant in a slowly varying magnetic field. As the field varies 



in time, the Larmor orbits expand or contract and the charged particles 
lose or gain transverse energy through change in velocity and angular 
momentum perpendicular to the magnetic field. This transfer of energy 
between particles and the field can be used for plasma heating, known 
as adiabatic compression. Since devices like TFTR include adiabatic 
compression in their design, transport codes used to predict TFTR's 
performance must contain the compression module. 

Alpha Particle Heating . An ignition reactor will be heated by alpha 
particles. Alpha particles are produced in the core of the reactor 
(50% of the a-energy is produced in 25% of the reactor radius). In 
order to determine the effect of a-particle heating, trajectories, 
confinement, loss cone, and slowing down mechanism for a-particles 
must be calculated. Determination of how a-particles diffuse and deposit 
their energy will affect temperature profiles of a reactor plasma. For 
example, if a-particles exhibit anomalous behavior due to instabilities, 
they will diffuse outward and will exhibit enhanced slowing down and . 
a flatter temperature profile will result. 

The modeling of classical alpha particle transport is in good shape. 

The effects of anomalous a-particle transport are not well understood. 
Theoretical studies in this area will receive increased emphasis as 
tokamak temperature increase. 

Plasma-Wall and Plasma-Limiter Interaction . Impurities produce radiation 
losses in plasma. Heavy impurities have led to hollow temperature 
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profiles due to large radiation losses. Temperature and thickness of 
the cold plasma layer near the wall are important in predicting impurity 
recycling rates and charge exchange processes. Particle and energy fluxes 
to the limiter and the walls depend among other things on the temperature 
distribution in the plasma column, impurity species, impurity density, 
density gradients, and microinstabilities. Impurity production and 
transport are an important part of calculating confinement times in 
fusion devices. At present, impurities are treated by prescribing a 
source for impurities based on estimates of impurity production. With 
higher temperatures and densities, heat loads on the walls and the limiter 
may lead to enhanced impurity production. This will have a major effect 
on the operation of a fusion device. The area of plasma-wall and plasma- 
limiter interaction has not been receiving sufficient attention in the 

past. Theoretical efforts exist and new studies are being initiated but 
more attention is necessary. 
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V. Recent Developments in Transport Codes 

The transport codes have been in a state of limbo for the past several 
years while the neutral injection tests were prepared. Concerning 1-D 
transport codes it was noted in 1975 that 


"Many models have been proposed to fit the data, and prompt 
the observation that these codes have reached maturity. 
Their adequacy can be tested against new experiments in new 
regimes, and further improvements may require a deeper 
treatment of the underlying physics." 


With respect to prospects for further development of transport calculations, 
it was noted that 


"...the future prospects for transport code development seem 
to be along the lines of multi-dimensional calculations in 
physical and velocity space. When considering the upgrading of 
geometry (two dimensional), beams (Fokker-Planck equation), 
impurities (many data files for all species), neutrals (neutron 
transport codes), and plasma transport (instability simulation), 
it can be seen that complexity can grow easily." 


In the interim period, development along the lines indicated above has 
occurred. There are now several 2-D transport codes, as well as two 
new efforts at ab initio calculation of transport rates (Grad-Nelson, 
Jardin-Hirshman) in 2-D geometry. There are also several Fokker-Planck 
transport codes and one code which does both 2-D and Fokker-Planck (to 
simulate TFTR injection/compression scenarios). Impurity radiation 
data files have been prepared for a substantial portion of the periodic 
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Table of the Elements, and impurity diffusion codes dealing with the 

simultaneous interaction of many charge states are now running. Monte 

Carlo and S treatments of neutral gas evolution have been developed as 
n 

transport code modules. Multi-dimensional calculations of the sawtooth 
(m-1) MHD instability and of the saturation of m > 1 instabilities 
to form helical islands have been reduced to prescription form for 
inclusion in 1-D codes. Detailed comparison with experiment is under- 
way, and will serve as a reference point for new results expected at 
high pp 

From the above, one can see that the code development has proceeded 
in the predicted direction. The codes have now been tested against 
new experiments in new regimes: the recent PLT neutral beam injection 
results have proven out predictions of the neoclassical ion confinement 
model long used in transport calculations. While overall energy confine- 
ment is regarded as "anomalous," there have been new developments in 
calculation and measurement pf high Z radiative effects and in the 
analysis of the role of large scale MHD islands on transport. Behavior 
previously seen as anomalous has subsequently been found to be dominated 
by radiation and/or MHD effects in some experiments. Since these processes 
can be modeled, there now exists a realistic prospect of resolving all 
anomalies. 

Future Directions 

i 

With the clear endorsement provided by recent experimental results, 
transport code workers now have the responsibility to create accurate 
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computational models with less regard for computational size or com- 
plexity. There is a new generation of high beta, non-circular, and 
divertor tokamaks which will begin operation soon, and the understandabl e 
features of low circular plasmas must be computed as accurately 
as possible to serve as a reference. (The plasma near the magnetic 
axis will be stable to ballooning modes, have large-aspect ratio, and 
be approximately circular). Some examples of such computational improve- 
ments include coupling large impurity diffusion calculations to transport 
codes, providing 3-D neutrals modules to calculate jet effects of gas 
puffing and of neutral beam sources on toroidal average neutral density. 

The coupling of MHD instability codes with transport calculations will 
also be useful in computing the effects of ergodicity when many modes are 
present and to calculate the effects of ballooning instabilities at high 
beta. Fully 2-D codes are essential in understanding the dynamic interaction 
between realistic transport-induced pressure and current profiles, and 
experimental shape control and MHD stability. This area of research, namely 
the use of 2-D transport codes in conjunction with MHD bal looning codes, 
will receive increased emphasis at each of the major laboratories. 

As the use of more physics and higher spatial and velocity dimensions 
increases in the future, codes will become increasingly specialized to 
particular problems or machines, e.g,, doublets, divertor tokamaks, 
etc. Through specialization, these codes will permit more efficient use 
of computational resources and will provide more detailed answers to 



16 


specific problems. On the other hand, development will continue on 
sophisticated and complex codes that will contain a large number of 
physics modules and will be used primarily to check the output of 
more specialized codes in order to validate the assumptions made. To 
take an example from the present use of this technique, 2-D codes which 
follow the detailed time-dependent motion of the magnetic flux surfaces 
are not necessary to adequately calculate loss of confinement due to 
collisions and microinstabil ities. Simple "fixed shape" approximations 
(even simple unphysical shapes such as ellipses) are more than adequate. 

The development of codes will be affected by the amount and quality 
of experimental data available to them. A more detailed discussion 
of this problem is given in the section on "Tokamak Data Reporting." 
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VI . Present Status of Physics Models 

There is a general consensus among workers in the field that some 
phenomena contributing to transport are well understood while others 
are not. Basically many areas dealing with atomic physics and neutrals 
are well in hand, but production and inward transport of recycled 
material is still poorly understood. Calculations for neutral beam 
deposition, alpha particle heating, plasma compression, and general gross 
2-D motion of the magnetic flux surfaces have reached advanced stages of 
development. 

The physical boundary region of the plasma is not understood quan- 
titatively. The partition of energy flow between wall and limiter depends 
on a detailed knowledge of the transport coefficients as well as non- 
coronal radiation from impurities. On the other hand, the source of 
the impurities depends on the partitioned energy flux and detailed 
knowledge of irreproducible surfaces.* Semi -quantitative understanding 
of wall reflux physics is important to reactor studies and long-pulse 
tokamak operation. For this a "systems approach" rather than a partial 
differential equation approach may be most appropriate. 

There is little understanding of the anomalous electron heat and plasma 
transport. This is due in part to the fact that experimental information 
from ohmically heated machines does not give a clear determination of 


*$ince it is not possible to exactly reproduce the condition of the walls 
and limiter from one shot to the next, we call these irreproducible 
surfaces. 
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the temperature dependence of the transport loss. Strong beam heating 
or a drastic radiation cooling of an ohmic discharge with a direct 
measurement of the conduction is required for this determination. Since 
temperatures from ohmically heated clean discharges are correlated with 
current and size for example, the absolute dependence on these 
variables cannot be determined without the knowledge of the temperature 
dependence. The simplest phenomenological description of anomalous 
transport is simply to assign a value to the anomalous conductivity 
proportional to 1/n. This gives Alcator density scaling. However, the 
value of the conductivity is not truly universal and a factor of 
3 to 5 variation from machine to machine has not been successfully 
correlated with any systematic dependence. It should be clear that such 
imprecision of experimental knowledge allows a number of crude theoretical 
models (with one or more adjustable coefficients) to claim consistency 
with a large part of the data. 

A major area of physics which has been studied for 2-D transport is 
the possible distortion of the plasma equilibrium by high- power neutral 
beam injection. In addition to heating and the formation of a high-energy 
tail in the ion-distribution function, high power neutral beam injection 
produces a toroidal current that affects the flux surfaces. The highly 
anisotropic spin-up of the plasma column also contributes to anisotropic 
pressure distribution (P„ f PjJ. The distortion of the flux surfaces 
due to pressure anisotropy can be handled with the generalized version 
of the Grad-Shafranov equation. None of the existing 2-D transport 
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codes can handle anisotropy or spin-up. The effect of anisotropic 
pressure should be comparable to the distortion of the surfaces from 
circular by the finite plasma pressure. The anisotropic pressure also 
implies a density variation along the flux surfaces of the order of the 
inverse aspect ratio. This density variation invalidates the usual 
argument for reducing the full 2-D transport problem to a 1-D flux 
surface averaged calculation. Due to great complexity of doing the 
problem correctly, the density variation effect on the transport may 
have to be ignored. 

A. 2-D Transport: Why and How 

A two-dimensional (2-D) configuration is defined to be one in which 
the flux function tp(x,z) (magnetic surfaces are defined by tp = constant) 
depends on both the radial and vertical cylindrical coordinates (x,z). 

The concept of 2-D becomes particularly important when the surfaces 
change shape in time. The increased interest in low aspect ratio, non- 
circular cross section finite beta tokamaks calls for the two-dimensional 
treatment of flux surfaces. Among the practical situations where geometry 
and finite beta can markedly alter transport rates and energy balance, 
are adiabatic compression, magnetic axis shift and its effects on neutral 
beam penetration, surface shaping to reduce MHD instabilities, pc 
divertors, rapid FCT heating, and low aspect ratio devices 

Since transport of ionized species (as opposed to neutrals 
lines is usually much faster than across field lines, tr« 
or even three dimensions can be reduced to a one-dimensional 
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dealing with transport across flux surfaces, augmented by a single 2-D 
equation describing the shape of flux surfaces. This approach has been 
called l^-D. Self-consistent I's-D transport calculations in which the 
transport processes (flows and surface-averaged transport rates) are 
computed self-consistently with the evolving surface shapes and magnetic 
field are necessary to obtain realistic pressure and magnetic field 
profiles to be used in magnetohydrodynamic stability calculations. 

The 1 description is valid for slow transport processes and excludes 
tearing modes, ballooning instabilities, and other fast processes that 
may remain inherently three-dimensional. 

There are at least two reasons for requiring two-dimensional rather 
than one-dimensional treatment of plasma transport. First, the 2-D 
(non-circular) geometry of the flux surfaces modifies the value of the 
plasma transport coefficients. The resistive transport coefficients 
for diffusion normal to magnetic surfaces involve surface averages of 
various magnetic field quantities which are sensitive to the variation 
of | B ] along J3. Thus, in high- e equil ibria, where \fy\ contours tend to 
align with the magnetic field, the neoclassical transport coefficients, 
as well as anomalous coefficients driven by magnetically trapped particles, 
could be substantially reduced from their usual 1-D values. For example, 
the vertical elongation of the plasma cross section increases neoclassical 
confinement time by decreasing the transport coefficients below the 
circular value by the square of the ellipticity. 
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Secondly, knowledge of the two-dimensional configuration of the plasma 
is necessary to apply boundary conditions and to compute the effects on 
the plasma of externally maintained heat and particle sources. Magnetic 
boundary conditions must be ultimately related to the currents in the 
poloidal field coils, whose spatial distribution is inherently two 
dimensional. The location of pressure contours, which approximately 
correspond to density contours, is important to determine heating and 
particle deposition profiles due to neutral beam injection. 

There are various degrees of complexity in the formulation of the lig-D 
transport problem. In the simplest 1VD technique the flux surface 
shapes are assumed to be known and fixed in time and area weighting of 
transport coefficients is used in standard 1-D transport equations. 

This method is often used for low beta plasma with fixed surfaces. 
Unfortunately, it is not known whether simple area weighting correctly 
reflects the change in transport rates due to surface shaping especially 
for the empirical jmode Is used to describe the electron dynamics. In 
general, for reasonable surface shapes and parameters confinement times 
calculated with this formulation will be within a factor of two of 
those given by a 1-D cylinder model. 

An intermediate technique goes beyond simply postulating the effects 
of geometry, but does not entail computing the complete transport flow 
problem. In this approach the surface-averaged fluxes are assumed to 
be known or postulated (from neoclassical theory or from a specific 



22 


model for anomalous transport) and they are used to calculate the change 
in plasma variables due to transport across surfaces of constant 
toroidal flux. As the plasma variables evolve, the 2 -D (Grad-Shafranov) 
equation is solved numerically to obtain the shape of flux surfaces. 

This new geometric information is passed to the transport equations. 


At present, the neoclassical transport coefficients needed to use this 
technique are not available in all collisionality regimes for arbitrary 
aspect ratio, arbitrary cross section or arbitrary 3. The near axis(es) 
region of low aspect ratio non-circular finite 3 tokamaks may be 
adequately described with present theory, since the aspect ratio there 
is large, the surfaces may be nearly circular and this region most 
stable to ballooning. However, volume averaged pis the parameter of 
interest, and this depends on global knowledge of the transport rates. 


Within the limitations of this model, however, there are experimental 
situations which it describes adequately. For example, the recent PIT 
neutral beam experiments employed a low n e , high T 0 target plasma to 
produce high proton temperatures. I n this case the ion and electron 
energy balances were decoupled (by low n 0 ) and the injection energy 
went predominantly to the ions (high T( , W . Adiabatic 

compression is another example where this approach is adequate, since 
the compression power supplies are chosen so that T . << 

Injection and compression coupled together are a 1 so'wen Elated , sincT^ 
the fast ion processes seem to be well described at 


present. 
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In the most sophisticated formulation of l^-D transport, the transport 
processes (flows and surface averaged transport rates) are computed self- 
consistently with the evolving surface shapes and magnetic field. As 
was mentioned before, the self-consistent 1^-D transport calculations are 
necessary to obtain realistic pressure and magnetic field profiles to 
be used in magnetohydrodynamic stability calculations. Once a stable 
initial equilibrium configuration is specified, the time-dependent trans- 
port equations determine how that equilibrium evolves in time. In 
particular, if the time-dependent equations lead to a steady state, this 
corresponds to a unique prescription for the equilibrium profiles. The 
importance of these resistive steady states in determining unique 
stationary profiles of the otherwise arbitrary pressure and toroidal 
magnetic field functions p(^) and F(i|>) has been previously emphasized 
in publications. 

Furthermore, a self-consistent 2-D calculation includes potentially 
important inductive contributions to the flux surface-averaged normal 
plasma flow arising from aJJ/at. The maintenance of the lowest order 
(in resistivity n ) pressure balance vp = ^ x J3 on the resistive time 
scale can lead to deformations of the plasma magnetic surfaces, with 
associated self-consistent inductive x B flows (JE^ - - aA/ at ) 
frequently greater than the resistive diffusion during the fast 
transient phase of the discharge. This inductive convection of the 
magnetic surfaces must be determined before the absolute particle flux 
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can be evaluated. The particle flux has a significant effect on the 
actual plasma confinement time and incorrect values for the particle 
flux will result in erroneous predicted confinement time. 


To date, this approach has been successful only for neoclassical transport 
in the col 1 isional (Pfirsch-Schluter) and collisionless (banana) regime* 
Transport coefficients in the more difficult intermediate (plateau) 
regime have not been obtained, except for large aspect ratio, elliptic 
cross section tokamaks, Some progress has been made in understanding 
shape effects on MHD and kinetic drift instabilities, which may cause 
anomalous transport. The self-consistent methods correctly calculate 
transport coefficients for the time-varying geometry, properly include 
convection and adiabatic changes and within the limitations of the 
neoclassical models treated, can handle any aspect ratio, surface shape 
(including separates), and beta with equal case. In the collisional 
case only 1-D transport equations of standard from plus a 2-D equation 
equivalent to the equilibrium equation are required. In the collisionless 

case integrals are also calculated on each flux surface to describe 
trapped particles. 


Even for the most sophisticated 1^-0 models presently in use the 

'« - ««. „ de . 

For many simulation problems only minor changes in transport will result 

from adopting a 1VD model and in these cases the additional complexity 

is probably not warranted. But for the class of problems described 
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above (and presumably for others) 1 3g-D models are already cost effective. 
Future problems such as transport from 3-D tearing or toroidal spin-up 
(due to momentum input from beam heating) may require a considerable 
increase in computer time. 

B. Neoclassical Transport Coefficients 
Although plasma exhibits a wide range of collective properties, most 
of this behavior takes place on a time scale shorter than the binary 
collision time. On a longer time scale, the so-called diffusion time 
scale, the plasma exhibits properties associated with collisional 
effects such as electrical conductivity, thermal conductivity, particle 
diffusion across or along magnetic field lines, etc. In order to 
calculate these macroscopic properties, transport coefficients must 
be obtained from the appropriate statistical description of the plasma. 
Transport coefficients are the quantities that linearly relate fluxes 
such as the energy flux, to forces such as the temperature gradient. 

When combined with the exact conservation laws for particles and energy 
and with Maxwell’s equations, the linear transport relations provide a 
closed set of equations that predict the time evolution of the plasma, 
given a proper set of initial and boundary conditions. 

The transport coefficients cannot be obtained by solving the so-called 
moment equations, since these are not a closed set of equations. For 
example, the first moment equation, called the momentum equation is ob- 
tained by taking an average over velocity space of particle velocities 
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weighted by the distribution function f(x,v,t). The equation for the 
first moment is of the form: 


m (n v) + vP - en {E + ~ v x B) * F 

'X 

where v is the mean velocity and P is "the stress tensor. 

Solution of this moment equation would give the macroscopic velocity of 
electrons at a given point in the plasma. But this equation contains 

'X ■> 

the stress tensor P and the friction force F. Closure relations for 

% 'X 

P and F in terms of the mean velocity are now available in all neo- 
classical regimes. However, viscosity and friction coefficients must 
be determined by solving for the distribution function. The fluxes 
can then be found by solving the fluid momentum balance equation given 
above. The distribution function for an ionized plasma obeys an integro- 

differential equation just like the Boltzmann equation for neutral 

-> ■+ 

particles except that it contains e E and v x B terms: 


Jf + V-Vf + S. (£ + I v x B) • * C(f) 

at nic 


where C is the collision operator that conserves particles, momentum 
and energy. Whereas In the Boltzmann equation, the collision term 
represents binary collisions, in a plasma, collisions are not pre- 
inary, because of the long range of the inverse square 
for Coulomb collisions. 
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The cumulative effect of small -angle collisions in a plasma is much 
more important than binary collisions. When a particle in a plasma 
is subject to simultaneous small-angle deflections, its progress in 
velocity space becomes a random walk superposed on the ordered motion 
due to macroscopic E and B fields. This yields an expression for 
the collision term of the Fokker-Pl ancle type: 

(lr) c - *7 |jV + 7 w: <VL) 5 c < f > 

where A.. and are "friction" and "diffusion" coefficients respectively. 

They are functions of x,v,t at which is to be calculated and 

1 c 

are also linear functionals of f (x,t). 

In transport theory applications, the Boltzmann equation is averaged 
over the rapid motion of a charged particle around a magnetic field line. 
The resulting equation, called the drift kinetic equation, has the form 

£ * (5„. - cm 

where v„ and £„ are the components of particle velocity and electric 

■+ 

field parallel to the magnetic field, and v^ is the drift velocity due 
to curvature o'f the magnetic field lines. 

The equation used to calculate ion thermal conductivity Is a drift- 
kinetic equation with the Fokker-Planck collision term that has been 
reduced to three independent variables: two velocity and one spatial 
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variable. One of the spatial variables is eliminated under the assump- 
tion of small poloidal gyroradius, the second spatial variable is 
ignorable because of toroidal symmetry. The third velocity variable 
is the angle of rotation around the magnetic field, which has been 
averaged out. 

Since the MHD time scale is much shorter than the diffusion time scale, 
in the calculation of transport coefficients the geometry is considered 
fixed and the transport coefficients are calculated for a given geometry. 

Present Status of Transport Coefficients 

The neoclassical ion thermal conductivity «. has been calculated only 
in certain limiting cases. Expressions for general values of collision- 
ality are known only in the large aspect ratio limit, while expressions 
for general aspect ratio and general geometry are presently available 
only in the collision-dominated and collisionless limits. In view of 
the importance of this transport coefficient in the interpretation of 
experiments on PLT and other large devices, for general geometry, aspect 
ratio, and coll isional ity should be calculated. 

As was mentioned above, the drift-kinetic equation is in three inde- 
pendent variables: two velocity variables and one spatial variable 
(assuming axi symmetry and small poloidal gyroradius). It would be 
possible to solve this equation, and calculate the required moments, using 
a large computer such as the CRAY. It would obviously take too much time 
to solve the kinetic equation during the running of the transport code 
for all of the significantly different flux surface configurations and 
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coll isional ities which occur. If all interesting geometries were 
specified by a small number of size and shape parameters, a table of 
values of k. could be prepared. With a minimum of three size and shape 
parameters , plus the coll isional ity, and assuming ten values of each 

4 

parameter, it would still be necessary to solve the kinetic equation 10 
times, and to make use of a table with this number of entries. 

C. Microinstabilities and Anomalous Transport 
One aspect of toroidal modeling that deserves attention is the influence 
of various microinstabilities on transport. The evidence from experiments 
points to the fact that the thermal transport of electrons is not 
governed by classical effects, Furthermore, there is difficulty in 
explaining the success of gas injection experiments on the b^sls of 
coll isional transport. In addition, high frequency wave scattering 
experiments reveal the presence of density fluctuations indicative of 
drift wave disturbances. On the basis of these facts, one is led to 
believe that these anomalous processes are the result of microinstabilities. 

The methods by which transport induced by microscopic fluctuation 
is modeled can be summarized as follows. It is assumed that the anomalous 
transport is the result of the appearance of some specific instability. 

A small amplitude non-linear or quasi linear analysis is performed to 
determine the relative size of the various transport coefficients. The 
magnitude of the fluctuations and the resulting transport are then 
related in some way to the degree of linear instability present. There 
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are areas of contention in each of the three steps just described. First, 
there is no universal agreement on which of the instabilities are 
important. Some of the modes thought to produce anomalous transport are 
the electron and ion drift modes, and drift-tearing modes. Secondly, 
the applicability of the methods by which one calculates the transport 
coefficients for a prescribed set of fluctuations has not been proved. 

For example, it has not been proved that electrostatic quasil inear 
theory is sufficient. If it is not, then the effect of fluctuations in 
the magnetic field is not known. Do the fluctuations lead to islands, 

i 

stochastic wandering of field lines, or some other effects? And 
thirdly, in the absence of some demonstrated saturation mechanism, it 
is not known how the level of fluctuations can be related to the degree 
of instability. One approach favors a weak' dependence on the linear 
growth rate in which case a precise determination of the growth rate 
is not needed. A second approach favors a rapid increase in fluctuation 
level as the profiles pass from a stable to unstable state, in which 
case, a precise determination of the growth rate is necessary. 


Evaluation of the various methods and models is difficult owing to 
the fact that there are few experiments in which these anomalous trans- 
port processes can be isolated from other losses and that it is difficult 
in general to measure the relevant unknown quantities. However, 
transport codes can be useful in testing the various models and’deter- 
mning their effect on confinement. Progress in this field would be 
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aided by experiments on smaller devices under more rigidly controlled 
circumstances and by particle simulation of anomalous processes. 

D. Treatment of Collisionless Ions in High Beta, Non-Circular 
Tokamaks 

Fast ions are characterized by their large deviations from a flux 
surface and by their consequent non-local nature. Thus, to properly 
treat the slowing down and plasma heating resulting from Coulomb 
collisions, the full orbit topology must be used. There are two 
ways that may be used to do this — a Monte Carlo particle following 
code, or a treatment of the fast particles in their constants of motion 

i 

space. 

Monte Carlo codes are relatively straightforward but are very time- 
consuming since the basic time step is a small fraction of a bounce time. 
Monte Carlo fast ion codes exist at Fontenay-aux~Roses , Princeton, and 
other places. To effectively use these codes, the slowing-down process 
must be artificially shortened and so far, less- than 1,000 particles 
have been used. If only moments of the fast ion distribution are required, 
this may give adequate statistics. But, there is some interest in measur- 

i 

ing the fast-ion slowing-down distribution function and comparing it to 
theory. This requires the use of many more particles than have been 
used until now. A unique advantage of the Monte Carlo approach is the 
ability to handle ripples in the magnetic field. However, this requires 
even shorter time steps and consequently longer computer times. A dis- 
advantage of the method is the lack of analytic results and hence it is 
very difficult to generate scaling laws or to predict trends. 
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The other approach to the problem is to express the Fokker-Planck 
equation in terms of the three constants of motion* which completely 
characterize the particle orbit in an axisymmetric tokamak. The 
bounce average of this equation describes how the constants of motion 
change due to collisions. This method is being explored at ORNL. 

The advantages of the constants of motion method are its slow time 
scale (much longer than a bounce time) and its analytic nature. Thus, 
there is hope that good approximations can be made and that scaling 
laws can be derived. The disadvantage of the approach is the difficulty 
of bounce-averaging the Fokker-Planck equations, a process which must 
be redone each time a plasma profile or equilibrium changes. 

Another approach is to use the "elements of the orbit" methods exploited 
in celestial mechanics. Here one computes a sample orbit using the 

constants of the motion. One computes the rate of change of the 

constants of motion and advances the orbit accordingly, thus avoiding 
the problem of bounce averaging the Fokker-Planck equation, or even 
solving the Fokker-Planck equation at all. 


*0ne such set of constants of motion would be: speed of the particle, 
at a point, maximum value of flux function along the particle orbit 
and v„ = p at the same point. 


v 



APPENDIX I 


The Baldur Code 

Baldur is a 1-D tokamak transport code that has been written primarily 
for TFTR. The code features two hydorgenic species (D and T), two 
impurities, two species of neutral gas, neutral injection heating, 
a-particle heating, flexible transport models, and slow adiabatic 
compressions. 

The code uses empirical rather than six-regime transport. It contains 
finite alpha slowing down plus loss cones for alpha particles. Particles 
from the beam are included in the density calculation. Effect of 1/3 
and 2/3 energy component and the compression of the real velocity 
distribution function with injection represent an improved treatment 
of neutral beams. 

The physics modules in the code can be divided into two groups pertaining 
to energy (Fig. 1) and particle (see Fig. 2) balance. They include 
electron and ion heat transport, ohmic heating, neutral beam heating, 
alpha-particle heating, slow adiabatic compression, impurity radiation, 
and neutral gas cooling. Modules dealing with particle balance 
include hydrogen and impurity transport (two hydroben and two impurity 
species), neutral gas transport (two hydrogen species), sources due 
to thermal! zed beam ions and sources due to thermal i zed alpha ions. 
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Numerical Features 

The code uses predictor-corrector method and first-order extrapolation 
for the integration of nonlinear, parabolic, initial value partial 
differential equations. Monte Carlo techniques are used for neutral 
gas transport, neutral beam deposition, and alpha-particle orbits. 

Simple Fokker-Planck schemes are used to calculate neutral beam heating 
and alpha-particle heating. 

Equations are finite-difference with variable spacing of the zones. 

Below is a brief description of individual physics modules contained 
in the code. 

Transport Models . The following models can be included by input: 
neoclassical, pseudo-classical, Bohm, microinstabilities, empirical 
K, D ^ T^ {x * -1, y * 0) ^ (a/r) X 

Normally the empirical model is used, adjusted to fit PLT and 
Alcator confinement data. 

X e * 10^/n e * cm 2 /sec 
D % 10 17 /n e cm 2 /sec 
K.. -v neoclassical 

An arbitrary mixture of the above-mentioned models can be specified. 

The code can also use the local drift mode and the trapped-electron 
mode to determine marginally stable profiles for different instabilities 
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for different coll i sional ity regimes- 

Fig. 3 gives a comparison between the experimental values of T as 

e 

a function of n in Alcator and the values calculated with two different 
values of X . Fig. 4 gives a comparison between experimental and 
calculated energy confinement time as a function on rT for Alcator. 

Two Hydrogenic Species . Separate transport equations are solved for 
each species. Separate neutral gas treatment is provided for each 
species, with charge exchange- between them. Thermal i zed beam ions 
are added to the plasma density. 

Two Impurity Species . The module includes coronal equilibrium radiation 
for 47 possible impurities (H to U with Z = 92). Separate transport 
equations for each impurity are solved in the course of the calculation. 
The module contains an optional simple non-coronal equilibrium. Fig. 5 
gives the power loss due to various impurities as a function of electron 
temperature . 

Neutral Gas . A Monte Carlo algorithm calculates transport of two neutral 
gas species. Neutral gas sources consist of recycling of plasma off 
the wall, gas puffing (edge effect), charge exchanged neutrals from 
the beam (volume effect) and recombination (volume effect). The 
assumption is usually made that gas puffing and recycling produce 40 eV 
neutrals. Fig. 6 illustrates particle tracking geometry used in neutral 
transport module. 
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Neutral Injection , The energy deposition profile is calculated with 
the aid of a Monte Carlo algorithm, A simple multigroup Fokker-Planck 
scheme is used to calculate the slowing down and pitch-angle scattering 
of beam ions. Slow adiabatic compression of beam ions includes pitch- 
angle effects. Loss of beam ions due to charge exchange is included 
but the recapture is not. The module contains no orbit effects and no 
beam-beam collisions. Fig. 7 shows calculation of tangential injection 
in TFTR. Fig. 8 shows a calculation of the time rise of T. and T 

1 G 

with beam injection. The figures show the result of Fokker-Planck 
and multigroup calculations. Fig. 9 shows electron heating in PLT 
before and after injection as a function of radius. 

Alpha-Particle Heating . A one-group Fokker-Planck scheme is used to 
calculate the energy deposition. Alpha-particles from beam-plasma and 
thermal reactions are included. Thermal i zed alphas are added to the 
plasma as an impurity. Loss cones are included for the first orbit. 

Loss cones and time-delayed heating are computed with Monte Carlo and 
Fokker-Planck schemes. Figs. 10 and 11 give the effect of 1^ and 
the aspect ratio respectively on prompt alpha-particle losses. Fig. 12 
gives density profiles for the creation of alpha-particles, the density 

of alpha-particles after losses have occurred, and the energy deposition 
by alpha-particles. 

Slow Adiabatic Compression . The compression is done at the end of each 
time step. Beam ions are compressed but alpha-particles are not. 

The limiter can be expanded during compression to study the gas blanket. 
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APPENDIX II: NYU IfrD TRANSPORT ALGORITHM 


Introduction 

The characteristic feature of the 1%~D family of algorithms is the 
separation of the solution process into alternate marching of profiles in 
!-D (current, temperature, flux, volume, etc.) and 2-D evaluation of 
contours and geometrical coefficients (surface inductance, surface 
averaged resistivity, etc.). The basic present limitations of the 
1%-D technique is to macroscopic, scalar pressure, MHD. Codes have been 
constructed and operated which include various combinations of finite e, 
general plasma cross section and aspect ratio, various topologies 
(Doublet, tearing, re versed-field mirror) including time-dependent 
transitions in topology due to external coil variation and plasma 
transport, with models including (classical) tensor resistivity and heat 
flow as well as the adiabatic limiting case. 

The three most important advantages of these methods over fully classical 
(including inertia and wave motion) and the very recently developed 2-D 
marching of the Grad-Hogan model are: 

1. Complex plasma geometry and topology, including dynamic changes in 
topology are easily managed; 

2. The code logic is such as to be easily adaptable to convert any 1-D 
(cylindrical) "Kitchen Sink" transport code to a fully 2-D basis. 

Extension of 1%-D techniques to closed line C 
completely 3-D toroidal transport seems nov 
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Density profiles for alpha births., alphas after losses 
and alpha energy deposition 



APPENDIX II: NYU 1%-D TRANSPORT ALGORITHM 


Introduction 

The characteristic feature of the 1^-D family of algorithms is the 
separation of the solution process into alternate marching of profiles in 
1-D (current, temperature, flux, volume, etc.) and 2-D evaluation of 
contours and geometrical coefficients (surface inductance, surface 
averaged resistivity, etc.). The basic present limitations of the 
l^-D technique is to macroscopic, scalar pressure, MHD. Codes have been 
constructed and operated which include various combinations of finite & 
general plasma cross section and aspect ratio, various topologies 
(Doublet, tearing, reversed-field mirror) including time-dependent 
transitions in topology due to external coil variation and plasma 
transport, with models including (classical) tensor resistivity and heat 
flow as well as the adiabatic limiting case. 

The three most important advantages of these methods over fully classical 
(including inertia and wave motion) and the very recently developed 2-D 
marching of the Grad-Hogan model are: 

1. Complex plasma geometry and topology, including dynamic changes in 
topology are easily managed; 

2. The code logic is such as to be easily adaptable to convert any 1-D 
(cylindrical) "Kitchen Sink 1 ' transport code to a fully 2-D basis. 

Extension of 1^~D techniques to closed line 3-D (EBT geometry) and 
completely 3-D toroidal transport seems now possible. There is, however, 
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no consensus among the participants of the workshop on the degree of 
difficulty in accomplishing this task. Strongly anisotropic problems 
(such as guiding center or drift models for a mirror machine with loss 
cone) do not seem amenable at present. The neoclassical (guiding center 
or drift but approximately scalar pressure ) problem seems to have 
reached a stage close to that of the okssical MHD Grad-Hogan model some 
five years ago. The theory including consistent geometry changes seems 
to be visible but algorithms and operating codes for a 1^-D or 2-D 
version need further technological improvement. 

The Basic 1%-D Algorithm 

As mentioned before, one part of the code is an equilibrium solver. 

In 2-D, given p(tf>) and B z = f ( ip) profiles, equilibrium is found by 
solving a nonlinear elliptic PDE, 

0) F(i (i), F = -dp/d* -fdf/d* 

One method of solution is to iterate using a Poisson solver, 

(2) Al Vl = F t* n (x.y)] = f n (x,y) 

This may or may not converge in cases where the solution is unique as 
well as when the solution is non-unique (i.e., after bifurcation). 

There is a large literature on more and more sophisticated procedures, 
but there is no universally dependable algorithm. In a case of multiple 
solutions, numerical stability (distinct from physical stability) 
determine which of the solutions may be found, 
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A second formulation of the equilibrium problem is 

(3) Ai|; = F(V) 

where the r.h.s, is a known function of V. The same iteration 

(4) A* n+1 = F(V n ), V n = v n (x, y} * * n (x,y) 

is usually more rapidly convergent than (1) but requires contours of 
and volumes to be evaluated. The numerical stability of solution 
algorithms for (3) is different from (1) and both are different from 
physical stability. 

A third equilibrium formulation* given adiabatic profiles v»U) and 
v { ip) defined by 

(5) P «* (V) ] Y 

f * v^' (V) 

leads (fory - 2 and v = 0) to the differential equation 

(6) ^ (<r) 2 - 2ur 5 r 

The solution will amost never converge on simple iteration 

A Vl = R n , 

(and this gives no opportunity to impose proper boundary conditions). 
A method which converges extremely well is to alternate between the 
Poission equation 
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(8) A Vl = R n (x ’ y) 

which is used to compute only the geometric contours ^ n+ i = const. 

(the ^ n+ i profile is discarded). Given the contours, the inductance 

(9) K n (V) = < |W| 2 > 

is calculated and the average pressure balance 

(10) (KJi ' ) ' = - ^(t') 2 - 2pf 

is solved from an ODE for ^(V), given [a priori] and given 

K(V) [temporarily, by iteration]. From ^(V) and the previous contours, 

V(x,y), one evaluates R n+ -|(x,y). 

Incidentally, the ODE allows specification of ^ at the center (V-Q) 
as well as the elliptic boundary condition at the plasma edge. 

In a non-adiabatic problem with transport, the averaged (1-D) transport 
equations take the place of Eq. (10), In 2-D, a representative transport 
equation 

+ u*V^ * hAiJj 

when averaged becomes 

+ Ity* * n (Kif> 1 ) 1 
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where 

9ii; 9 a 

3t = It ^ (x,y, t) , t|» t = if>(V,t) 

and 

U = <u-w> = j u . ds 

If an adiabatic variable (e.g., mass) is taken as an independent variable 
+ Uty’ » dty/dt « \|>(M,t) 

the velocity U (the plasma diffusion velocity) is eliminated and the entire 
set of transport equations becomes a set of 1-D transport equations in M 
and t, with various geometrical coefficients similar to K. 

There is always one fewer averaged transport equation than originally 
(since one of’ the “diffusing" variables is taken as the independent 1-D 
variable). The count is made up by adding the average pressure balance 
(10) [in an appropriate independent variable, not necessarily V]. 

The skeleton algorithm is: 

1. From initial equilibrium calculate (contour averaged) coefficients for 
1-D averaged transport equations. These are marched in 1-D, together 
with the average pressure balance for many diffusion time steps for 

one geometry time step. 
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2. Using the diffused 1-D profiles at the advanced geometry time step, 
a r.h.s. R(x,y) is computed for the Poisson equation. [Here one has 
a choice, involving simplicity, accuracy, stability, etc. of using 
the forms (1) or (3) or (6) since aJJ_ profiles are available from 
the result of diffusion]. 

There are two iteration loops, the inner one to solve for a geometry 
F(ij>) or F( V) or p(<j>), and the outer one which repeats the diffusion 
over a geometry time step until all geometrical profiles such as K(V) 
converge. 

One can take K(V,t) to be piecewise constant in t, or make it linear 
in t between geometry steps (for insertion in the transport equations), 
or more sophisticated (e.g., predictor-corrector). 

There are at least three (and preferably four) independent meshes 

1. 2-D - (x,y) 

2. 1-D - contours in (x,y) to evaluate averages 

3. 1-D - diffusion 

4. 1-D - interpolation among other meshes. 

iicn tw/n semi -independent convergence criteria and cq 

’Jter loops. The optimization problem is formidable 
oal. For example, in studying the small resistivity 
c, there is a very sharp moving singularity in K(V,t) 
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at the separatrix; this requires careful non-uniform contour placement. 

To optimize computer time, more sophisticated extrapolations of K(V,t) 
in t must be used, varying the values of and eq, the relative sizes 
of 2-D contour, and diffusion meshes, must be performed to reduce the 
running time. In one problem at NYU dealing with transition from a Belt 
Pinch to a fully formed Doublet with practical parameters, proper use 
of contours and diffusion meshes made it possible to run the code with 
only five updates of the geometry. 

In a transition from circular to D-shaped (or change induced by increasing 
beta) it should rarely be necessary to use more than two or three 
geometry calculations; but this will have to be confirmed in different 
parameter ranges. 

More details of the equations and algorithms are found in the following 
publ ications : 


1. Grad, H., Hu, P. N. and Stevens, D. C, "Adiabatic Evolution of 
Plasma Equilibrium," Proc. Nat. Acad. Sci.» USA, 72, pp. 3789-3793 
(1975). 

2. Grad, H., Hu, P. N., Stevens, D. C. and Turkel, E., "Diffusive 
Transition from Belt Pinch to Doublet," 2nd European Conf. on 
Computational Physics, Computing in Plasma Physics and Astrophysics," 
Garching, Germany, April 27-30, 1976. 

3. Grad, H., Hu, P. N., Stevens, D. C. and Turkel, E. "Classical 

Plasma Diffusion," IAEA 6th Conf. on Plasma Physics and CTR, Berchtesgaden, 
Germany, Oct. 6-13, 1976. 

4. Grad, H., Hu, P. N. "Classical Diffusion: Theory and Simulation 
Codes," Workshop of High 3 Plasma, Varenna, Italy, Sept. 1977. Also 
Rept. MF-91 , DIMS, NYU, March 1978. 

5. Oak Ridge National Laboratory, Theory Department Memo #75/9. 
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A proof of existence (not convergence of the numerical algorithm) has 
been given for a linearization of the GDE fG. Vigfussen, Thesis, 

% 

NYU, 1 977] which is the basis of the 2-D marching algorithm employed by 
Jardine (PPPL) for the Grad-Hogan model (these linearized versions are 
even more useful for analytical work than numerical). 



APPENDIX III: A 2~P RESISTIVE TQKAMAK AND DOUBLET TRANSPORT CQD £ 


Purpose : 

1. To construct an efficient ”1 1/2 D" transport code in 
the Grad-Hogan formulation. 

2* To investigate the range of validity of Pf irsch-Schluter 
diffusion . 

J > . To investigate the coupling between geometry and skin 
effect, geometry and plasma diffusion, skin effect and 
plasma diffusion. 

4 . To construct an efficient MHD code with real geometry 
as a foundation for next generation "Kitchen Sink" 
Tokamak and Doublet simulation. 

To compare resistive with adiabatic transition from 
Doublet to Belt Pinch. 

6. To investigate current profile evolution resulting from 
geometry changes. 

7 * To investigate plasma compression on the skin time scale. 


XI Description 


1. Low f3, scalar ■ pressure, classical 
resistivity, MHD. 

2. 2D, plasma separated from rectangu- 
lar box by vacuum region, with coils. 

3 . Toroidal and poloidal external 
field individually varied. 

4. Primarily on resistive skin time 
scale (slower plasma diffusion 
superposed) . 

5» Tokamak, Belt Pinch, Doublet, and time dependent transitions. 

6. Ability to handle low resistivity (adiabatic transition) 
with accompanying large, localized current layers (at 
plasma edge and at separatrix). 

7* Plasma edge BC with skin effect and skinless. 

8. Specialized to y = 2. 







Ill* Equations 


Low p formulation [Nice, 1970] 


4* u*7^/ = rfii/ - c(t) 


A^ - [J 


= UjK 

U Q « a(t)V , a 



df 
o 

dt 


U is the plasma compression induced by changing zero 
o 

order toroidal vacuum field, f (t). 

The net plasma diffusion , = 
order on the skin time scale and can be commuted (separately) 


by 

the formulas 

V 

f 

V 


«i e v K *'*t 

~ V av - 

V(J^‘)^dV [Grad-Hogan] 

or 


V 

4 

0 


u p “ - nKp' - 

CK^' 

[Pfirsch-Schluter} 


The complete formulation as used for computation is 

+ V’ m n C * ) * -c 


jf u«ds (= <U‘VV>) is higher 


A^ « P 




F = ( 1 ) 1 



IV . Boundary Conditions 


1) A poloidal "boundary condition f = ^(t) [usually constant in 
space, but not necessarily] is applied on the external 
rectangular box, 

2) Coils (in this code, symmetrically placed) carry . varying 
currents I (t) which also influence the poloidal field. 

3) The toroidal vacuum field f (t) can be varied thru the 
convection term U . If f is constant, the plasma volume 

V is constant [otherwise V f ~ const.]. 

P Jr U 

4 ) In this scaling, poloidal field pressure is not strong 
enough to compress the plasma; changing if (or I ) will, 
however, induce a skin current unless a special skinless BC 
(taken from the adiabatic pr'oblem) is used. in most cases, 
one of the following BC's has been used 

a) Ip ~ given at box 

b) - const, at V - Vp 

c) = const, at V » 

BC (b) (constant total plasma current) gives a skin current 
layer; BC (q) can be shown to coordinate (higher order distinct 
from f ) toroidal and poloidal knobs so that there is no skin 
boundary layer. 

5) At V - 0 a conventional heat equation natural BC is imposed 
on 

6) Treatment of Separatrix is described separately. 



V. Basic Algorithm 


For a given K(V, t), rp is inarched by Solving a conventional 
ID diffusion equation in 0, < V < V with conventional BC at 

V *= 0 and choice (usually of (b) - skin or (c) - skinless) 
at V . 

Xn an early version of this code [results presented at 
Madison, 1976], the separatrix was ignored and the two symmetric 
islands were treated as described in Code 1, taking 

V - ■+ V 2 ~ 2V^, as independent ID variable and assuming that if/ 

and Kif / 1 are continuous across the Sep. (essentially by ignoring 
its location or existence). A later, more accurate version is 
described in Sec. VI. 

K(V,t i ) is calculated in 2D at intervals by taking 
J ( V ’ t^) J(x,y, t. ) from the diffusion equation, inverting 

£ “ J(x,y) and computing a new K(V, t^) from the contours of ip ^ . 
Between t^ and K is assumed to be linear in t. The outer 

^°°P j assigning that all quantities are knovm at t^, diffuses 
to recomputes K i+1 , diffuses again from to with 

the new K until convergence, 

The inner loop is an iteration at a given time t^, given 
J(V), in which ip ^ and the contours are iterated, kip ® J[V(x,y)], 


V(x,y ) rp(x,y ) [Poisson] 

^(x,y) *^V(x,y) [volume within contour] 



VI. Separatrix Singularity and Normalized Variables 

At the separatrix, K is unbounded as log |V-V S |; this is 

nontrivial since K'~ V|V-V B | occurs in the inhomogeneous term, 

j = (W) 1 , of the Poisson equation. The actual singularity 

can be evaluated analytically and is not quite so bad. For a 

moving singularity, K - log, in a diffusion equation 

^ thg current density is bounded but has a cusp almost 

t 

at the separatrix (preceded by a smooth dip, it T) is sufficiently 
small). The height of the cusp is unbounded as rj becomes small. 

The following normalized independent variables are taken, 
inside and outside V sep 

V = VA S 0 < V < 1 

V - (v-v g )/(v p -y s ) , 0 < V < 1 

In (V, t) the diffusion equation has a convection term, (VVV ♦ 

s s 

Both V g (t) and K(V, t) are obtained from the geometry at inter- 
vals t^, K by linear interpolation in t, and V (t) by a cubic 
spline (so that V g (t} is continuous at t^). 

A moving log is fitted to K near V « V g (t) by taking a pair 
of contours close to V on each side. K * a. +b. log |V~V | io 
to each side (i = in, out); is replaced by |( b x+ b 2 ) 
to the last contour on each side. The log fit 
ch difference with a very fine 2-D mesh, but It 
s accuracy with a coarse mesh (which cannot see 



VII. Meshes and Interpolation 


A fine mesh is used separately on each side of the Sep. 

for ID diffusion. K, evaluated on contours (a fixed number 

outside, varying with V D inside), is transferred to the diffusion 

s 

mesh by cubic interpolation. The fine diffusion mesh is also 
used for transfer of J(V) to the 2D mesh by linear interpolation. 
The contour spacing has several alternatives: 

1) number of contours and distance of first contour from Sep. 
give a geometric spacing (constant ratio of adjacent spaces), 

2) number of uniformly spaced contours specified, 

3) optional subdivision of 11 last” sibdivision (near V = 0 and 

near V = V ) . 

P 

Although possibly unreasonable at first glance, the combina- 
tion of ID and 2D meshes allows very fine resolution with 
moderate 2D mesh and contours (say 32x64 and 10 outside contours) 
in which a very sharp current profile of total width less than one 
2-D mesh is accurately resolved, and an accurate cusp is drawn. 

A number of algorithms for further improving the accuracy of 
description of small islands has not yet been fully implemented. 



The equilibrium solver is very sensitive where the current 
peak at the Sep. is very large (a small change in V g makes 
large changes in J(x,y)]. The outer loop is found to have 
properties similar to an asymptotic series when the current peak 
is large (small .rj and small island); viz. there is an optimum 
outer loop time step - k; if too large or too small, many 

iterations are required. 

In a transition from simple geometry to islatcd, the cock 

first notices a change in topology in the 2D Poisson inversion. 

If the island is too small, it is ignored (no contours are close 

and K(0) * 0 is used together with all finite, simple contours. 

When the island is large enough to fit three internal 

contours (plus V - 0 and the Sep. ) after convergence of the 

geometry at t^j K is computed on two normalized meshes and 

transferred to the diffusion mesh before diffusing to t At 

i+1 

the same time, $ is recomputed over the island region usinc 
(K^' ) 1 = J and the new K. 



VIII. 

1) 


2 ) 

Not yet 

a) 

b) 


Auxiliary Algorithms 

The outer loop is almost always oscillatory (this can 
be qualitatively justified theoretically). There is 
therefore a great advantage in convergence to back- 
average the bare result of iteration with the previous 
iterate. A simple exercise in control theory adjusts 
the back average coefficient in terms of the observed 
convergence rate. 

Back averaging is used in both outer and inner loops; 
the former on the ID array , J(v), and the latter on the 
2D array, $ . 

implemented. 

Refining 2D mesh in Box which surrounds island. 

Using analytic formulas for K(v) in small island, 

extrapolating back to precise time of islation with 

3 /2 

analytic growth ~ t 7 . 



IX. Results 


1. Sample results including a demonstration of enormous 
deviations from Pf irsch-Schluter were reported at the Sherwood 
Theory Meeting, 1975, and published in the Second European 
Conference on Computational Physics, Garching, 1976. 

2. Very accurate current response to deformation, 

including a current spike precursor before islation, transient 

skin effect signature at the center of the plasma from 

2 

oscillatory external coils, calculation of heating at 
Pi asma center from oscillating separatrix current, and accurate 
reversible (half sine wave) operation at small r) were reported 
at Berchtesgaden in 1976. 



(c| rally Dcvcloptfd island* 
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APPENDIX IV: TOKAMAK DATA REPORTING 


We feel that it is essential that there be an improvement in the 
methods of data reporting. We recognize that all the important 
experimental data are not always available on every series of discharges, 
Our major concern is that the reduced data is often not available in 
an organized, standardized, and easily accessible form for widespread 
interlaboratory use in transport codes. There is an apparent lack of 
understanding in many cases as to what information is generally useful 
for this purpose. 

In published experimental papers, the reduced data is often presented 
only in the context of analysis, interpreta tion, and experimental proof 
of a hypothesis. Thus only that part of the data thought to be relevant 
to the problem at hand gets reported. Information correlated on a 
shot-to-shot basis is not available in published papers. For example, 
it is of little use for transport code work to show "typical" profiles, 
time traces, or values of temperature uncorrelated with density, current, 
field, or working gas, etc. In many cases (e,g., scaling analysis), it 
is sufficient to report the "steady" parameters. Sometimes the ful 1 
dynamic behavior is of interest (e.g., gas >■ 

curves, information on whether or not the 

We feel that understanding of results o 

generation of large devices will be dramatically improved if the large 
laboratory experiment (in particular, PLT, Alcator, ISX, and Dill) have 
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scientific "secretaries ," as in the case of TFR, to collect and correlate 
the available reduced data in a standard form. What data should be 
reduced and in what form can only be decided by the experimentalists 
and interested theorists within the laboratory. This is already being 
done to some extent in each laboratory. The reduced data should be 
placed on the CTR network computer and made available to all interested 
workers. The time lag between data reduction and placement in the 
system should be sufficient to: (1) allow the data to "age" (i.e., 
until the experimentalists feel that it is reasonably correct), and 
(2) to allow some priority rights of the experimentalists to interpret 
the experiments "first." However, under the present conditions the 
time for data to become public is long compared to the time between 
major experiments. 

Below we give some examples of data often available but not routinely 

i 

reported in a correlated fashion. 

1. Measurement of the radiation power loss density versus radius. Avail- 
ability of this data is essential to understanding transport. If 
this quantity is significant relative to the total local power balance 
of the electrons, no treatment of the electron power balance can be 
given. It is insufficient to state the total power radiated; the 
discharge can be transport-dominated even when 100% of the power ' 
goes to radiation provided the central region has little radiation 
loss. Radiation (from high-Z impurities in particular) is uncorrelated 
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with any other parameter and because of its dependence on the 
surface source, it is a priori unpredictable in practice. 

2. The neutral density at the center. In conjunction with standard 
charge exchange codes this data allows some inference on the 
plasma transport coefficient. The plasma transport relative to 
electron heat transport is an important signature of all transport 
theories. The neutral density at the center can be obtained from 
the absolute calibration of the charge exchange diagnostic which 
measures the ion temperatures. 

In addition to mating reduced working data available for transport 
code work, we feel that there is a need for more confinement scaling 
experiments on lokamaks. By scaling we simply mean a series of dis- 
charges under comparable machine conditions in which only one variable 
is changed. Tin* reporting of "typical" discharges (frequently it 
means discharges having the maximum achievable confinement time) is 
of limited usefulness in transport codes since transport models 
will Invariably have at least one adjustable parameter. As the tokamaks 
become cleaner (or if the radiation can be subtracted out as described 
above) and the shots become more reproducible* such experiments will 
take on added significance. Since we are limited to factors of 2 
and 3 In many of the experimental parameters, fairly careful inter- 
pretation and analysis with transport codes is needed to correlate 
the data and account for profile changes due to radiation and MHD effects 
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A community-wide review of this question is needed. A standards 
committee is needed which would establish reasonable guidelines for the 
dissemination of that data which never sees publication in finished 
form, but which is vital to detailed analysis of confinement processes. 
The predominant membership of this committee should be experimental. 




